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Abstract 

^vq New versions and extensions of Benson's outer approximation algorithm for solving 

linear vector optimization problems are presented. Primal and dual variants are provided 
' . 1 in which only one scalar linear program has to be solved in each iteration rather than 

two or three as in previous versions. Extensions are given to problems with arbitrary 
pointed solid polyhedral ordering cones. Recent set-valued duality theories turned out 
(-h to be crucial for the developments. The need to compute the values of a set-valued risk 

measure for multi-variate random variables was a driving force for this work. Numerical 
examples are provided, one of them even involving a new set-valued risk measure for 
multivariate positions. 
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1 Introduction 

<N 

Set-valued approaches to vector optimization are promising in theory and applications. A 
1 duality theory in this framework is important for algorithms, and the dual problems can be 

interpreted in certain applications, see e.g. [13 [13 I2D (221 H21 HH ESI [HI El E2J- Benson's 
outer approximation algorithm is a fundamental tool for solving linear (and also convex) 
vector optimization problems [H] 01 [H EU [32]. It is also important for solving set-valued 
problems |27j . Recent applications of linear vector optimization concern financial markets 
with frictions (transaction costs). For such applications, one obtains optimization problems 
which are genuinely set- valued. 

In this article, we introduce a primal and a dual algorithm of Benson type where only 
one LP has to be solved in each iteration step. In contrast, previous versions [3[ I4"[ 18] l3T| 132] 
require at least two different LPs in each step. As the main effort of Benson type algorithms in 
typical applications is caused by the LPs, the computational time can be reduced considerably 
by the new algorithms. Another advantage is that all LPs have a very similar structure and 
therefore the impact of warm starts can be improved. A further benefit is an improvement 
of the error estimation given in |31} [32], i.e., in approximate variants of the algorithms: The 



same approximation error can be achieved with fewer iteration steps (compare Remark 4.10 



and Example |6.1| below). For both the primal and dual algorithm two variants ('break' and 



'no break') are presented and compared (compare Example 6.1). Another novelty of this 



article is that linear vector optimization problems with arbitrary polyhedral solid pointed 
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ordering cones are treated, whereas in all other references [31 IH [U EJ EH EU [32] only the 
special case of the usual ordering cone R 1 ^ is considered. This feature will be exploited in 
applications involving set-valued risk measures for multi-variate random variables in markets 
with transaction costs. In such situations, ordering cones are usually different from M.^_ and 
generated by a large number of directions. A short introduction into this topic and several 
(numerical) examples are given. Examples 5.2 and 6.6 involve a new type of a set-valued risk 
measure which we baptized the 'relaxed' worst case risk measure. 

This article is organized as follows. In Section 2 we provide some basic notations and 
results. All of the remaining three sections start with short introductions. Section 3 contains 
an overview on the set-valued approach to linear vector optimization and related duality 
results where, in contrast to most of the literature, we allow ordering cones more general 
than In Sections 4 we introduce the new variants of Benson's algorithm. We also give 
a detailed description of the two-phase-method to treat unbounded problems. Section 5 
provides an introduction to applications involving set- valued risk measures, and in Section 6 
several numerical examples are reported. 



2 Preliminaries 

Let A C R 9 . We denote by clA, int^4, bdA the closure, interior and boundary of A, re- 
spectively. The set A is called solid if its interior is non-empty. A convex polyhedron or a 
polyhedral convex set A in M 9 is defined to be the intersection of finitely many half spaces, 
that is 

r 

A = P| {y G R«| (zYy > lt } 
i=i 

where z , . . . ,z r G M 9 \ {0} and 71, . . . , 7 r G M. As polyhedra considered in this article are 
always convex, we will not mention convexity explicitly. Every non-empty polyhedron ACM 9 
can be expressed by means of a finite number of points y , . . . ,y s G M 9 (s G N \ {0}) and 
directions k 1 , . . . , k l £ R q \ {0} (t G N) through 

S S t S N 

M=l j=l 1=1 ' 

where k G M 9 \ {0} is called a direction of A if A + {A • k} C A for all A > 0. This can be also 
written as 

A = conv {y 1 , ... ,y s } + cone {k 1 , ... ,k f } . (1) 

Note that we set cone0 = {0}. The polyhedron A is bounded if and only if the cone-part in 
the above formula is {0}. The vectors {y 1 , . . . , y s } and directions {k , . . . , are called the 
generators of A. The set A^ := cone {/c 1 , . . . , A;*} is the recession cone of A. A finite set of 
half spaces defining a polyhedron A is called H-representation (or inequality representation) 
of A, whereas a finite set of points and directions defining A is called V-representation (or 
generator representation) of A. A bounded polyhedron is called a polytope. A convex subset 
F of a convex set A is called a face of A if (y, y G A A AG (0, 1) A Xy + (1 — X)y G F) implies 
y,y G F. A set F is a proper (i.e. / F / A) face of a polyhedron A if and only if there is a 
supporting hyperplane H to A with i* 1 = .ff D A. The proper (r — l)-dimensional faces of an 
r-dimensional polyhedral set A are called facets of A. A point y G A is called a vertex of A if 
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{y} is a face of A. If k G \ {0} belongs to a half-line face of a polyhedral set A, then A; is 
called extreme direction of A 

A polyhedral convex cone C C M q is called pointed if it contains no lines. Of course, a solid 
and pointed convex cone is non-trivial, that is, {0} C C C 1', A non-trivial convex pointed 
cone CCf defines a partial ordering <c on M. q by y 1 <c y 2 if and only if y 1 —y 2 G C. If C = 
I£+ := {2/ €E M. 9 ! j/i > 0, . . . , y q > 0}, then the component-wise ordering < r <j_ is abbreviated to 
<. A point y G R q is called C-minimal in A if y G A and ({y} — C\ {0}) n ^4 = 0. The set of 
C-minimal elements of A is denoted by Mine A. If C is additionally solid (but not necessarily 
pointed), a point y G R 9 is called weakly C-minimal in A if y G A and ({y} — int C) n j4 = 0. 
Likewise, by replacing C by —C, C -maximal and weakly C -maximal points in a set i C K' 
are introduced and we write Maxc A for the set of C-maximal elements of A. The dual cone 
of a cone C C R« is the set C+ := {y* G R 9 | Vy e C : {y*) T y > 0}. The £-th canonical unit 
vector in R 9 is denoted by e\ 

3 Linear vector optimization 

In this section we outline the set-valued approach to linear vector optimization including 
duality theory and establish a more general setting where arbitrary polyhedral ordering cones 
C QM q rather than C = M.^_ are supposed. A comprehensive exposition for the case of the 
ordering cone C = can be found in [25] . The origin of this approach is discussed in \25\ 
Section 4.8]. A related duality theory and an overview on other approaches to duality for 
linear vector optimization problems can be found in a recent paper by Luc |28j . 

3.1 Problem setting and solution concepts 

The solution concepts defined in this section are based on the idea that in vector optimization 
(in contrast to scalar optimization), minimality and infimum attainment are no longer equiv- 
alent concepts. In order to have an appropriate complete lattice where an infimum is defined 
and exists, a set-valued reformulation of the vector optimization problem is necessary. Here 
we just introduce the solution concepts that result from these ideas. We motivate these con- 
cepts from an application oriented viewpoint only. More details and a theoretical motivation 
can be found in [2B E51 E! • 

Let matrices B G M mxn , P £ M qxn , a vector b E M. m and a solid pointed polyhedral cone 
C C M q be given. The following linear vector optimization problem is considered: 

minimize P : M n —> M. q with respect to <c subject to Bx > b. (P) 

Define 

S = {x G R n I Bx > b} and S h = {x G R n | Bx > 0} . 

Of course, we have S h = Sqo, that is, the non-zero points in S h are exactly the directions of 
S. A point x G S is said to be a minimizer for ([P]) if there is no x G S with Px <c Px and 
Px 7^ Px, that is, Px is C-minimal in P[S] := {Px\ x G S}. The set of minimizers of ([P]) is 
denoted by Ming. 

A direction k G M n \ {0} of S is called a minimizer for ([P]) if the corresponding point 
k G S h \ {0} is a minimizer of the homogeneous problem 

minimize P : M n — > M. q with respect to <c subject to Bx > 0. (P h ) 
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Definition 3.1. A set S C S is called a set of feasible points for (|P]) and, whenever 5^0, 
a set S h C S h \ {0} is called a set of feasible directions for Q. 

A pair of sets (§, S h ) is called a finite infimizer for ([P]) if S is a non-empty finite set of 
feasible points for Q, S h is a (not necessarily non-empty) finite set of feasible directions for 
0, and 

conv P[S] + cone P[S h ] + C = P[S] + C. (2) 

An infimizer can be understood as a feasible element for a set- valued extension of problem 
( [P"] ) (lattice extension) where the infimum (which is well defined for this lattice extension) is 
attained, i.e., condition ^ stands for the infimum attainment. The set V := P[S] + C is called 
upper image of (|P]). Clearly, if (S, S h ) is a finite infimizer and ({0} ,Y) is a V-representation 
of the cone C, then (P[S], P[S h ] U Y) is a V-representation of the upper image V. 

The following solution concept is based on a combination of minimality and infimum 
attainment. 

Definition 3.2. A finite infimizer (S, S h ) of (|P]) is called a solution to ([P]) if S and S 

consist of only minimizers. 

In practice, the upper image V is one of the most important information for a decision 
maker. This is due to the fact that in typical applications the dimension n of the decision 
space is considerably larger than the dimension q of the outcome (or criteria) space. The 
problem to calculate all the minimizers is usually not tractable. Moreover, the overwhelming 
set of all minimizers is in general not suitable to support a decision. It is more natural and 
easier in practice to compare the criteria rather than decisions. Furthermore, also in scalar 
programming it is often not necessary to know all optimal solution. A solution as introduced 
above can be seen as an outcome set based concept which provides the information to describe 
the upper image V . 

3.2 Duality 

If vector optimization is considered in a set-valued framework, it is very natural to consider a 
dual problem with a hyperplane- valued objective function. First, we note that this is also the 
case in scalar optimization as in R a point and a hyperplane are the same thing. Secondly, 
we have in mind the well-known dual description of a convex set by hyperplanes. The upper 
image V of a linear vector optimization problem (|P]) is a convex polyhedron which can be 
interpreted as an infimum of the lattice extension of problem ([P]), see [25J. It is therefore 
natural to ask for a dual description of this convex set which is obtained as the supremum 
of a suitable dual problem. As a third argument, we want to mention that there is a lack of 
applications of the classical approaches to duality theory in vector optimization. For instance, 
Ehrgott |6j pointed out that " dual algorithms could not be developed because of the absence 
of a duality theory for MOLP that could be algorithmically exploited." 

The idea of geometric duality [19] is to transform the hyperplane-valued dual problem 
into a vector optimization problem. This idea is taken from the theory of convex polytopes, 
where an H-representation of a polytope A defines a V-representation of a dual polytope. For 
instance, if A is solid and contains zero in its interior, an H-representation of the form 

A = {xeR q \ Bx< (1,...,1) T } 

exists. The row vectors of the matrix B yield a V-representation of the polar set A° := {y* 6 
K*| Vy G A : y* T y < 1} of A, which is a dual polytope to the polytope A. The duality relation 
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between A and the dual polytope A° is given by an inclusion reversing one-to-one map $ 
between the set of all faces of A and the set of all faces of A°. 

A similar duality map can be used to transform a hyper plane- valued optimization prob- 
lem into a vector-valued problem, which is called the geometric dual problem. We assume 
throughout that there exists a vector 

c G int C such that c q = 1, (3) 

and we fix such a vector c. As C was assumed to be a solid cone, there always exists some 
c 6 int C such that either c q = 1 or c q = — 1 . In the latter case we can consider the problem 
where C and P are replaced by — C and — P, which is equivalent to (P) and fulfills Q. 
Consider the dual problem 

maximize D* : R m x R q -> R q with respect to < K over T, (D*) 

with (linear) objective function 

D* :R m x R q ^R q , D*{u,w) : = (w x , w q -t, b T uf , 

ordering cone 

K:=IR + -(0,0,...,0,1) T , 

and feasible set 

T := {(u,w) € R m xR q \u> 0, 5 T -u = P T w, c T w = 1, y T w > 0} , 

where Y is the matrix whose columns are generators of the ordering cone C. 

A point (n, w) G T is said to be a maximizer for (D*) if there is no (u,w) G T with 
D*(u,w) >k D*(u,w) and D*(u,w) / D*(u,w). The set of maximizers of (D*) is denoted 
by Max (TEFl) . 



Definition 3.3 ( |16l I20 | 125 1 1 1 7 j ) . A non-empty finite set T of points being feasible for (D* 
is called a finite supremizer of ( D* ) i/ 

conv D* [T] — K = D* [T] - K. (4) 



The set T>* = D*\T\ — K is called lower image of (D* ). Similarly to above, a finite suprem- 
izer yields a V- representation of T>*. Condition Q can be interpreted as the attainment of 
the supremum in a suitable complete lattice, see e.g. [25]. The combination of maximality 
and supremum attainment leads to a solution. 



Definition 3.4. A finite supremizer T of (D*) is called a solution to (D*) if it consists of 
only maximizers. 

Note that, in contrast to ([2]), we do not need directions in Q. This is due to the simplicity 
of the cone K in contrast to C. 

A duality mapping ^ for the vector optimization problems ([P]) and (D* ) is now introduced. 
The bi-affine function 



9-1 / 9-1 \ 

; ■ r:> r. • - r.. <p(y, y*) := ^ yiy* +y q ll — ^ c iy * - .//,, 

»=i V i=i J 
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is used to define the following two injective hyperplane- valued maps 

H : R q , H(y*) := {y G R q \ v?(y, y*) = 0} , 

H* : R q ^ R q , H*(y) := {y* G R q \ <p(y, y*) = 0} . 
The mapping H yields the duality map 

V:2 Rq ^2 Rq , *(F*):= f] H{y*)nV. 

y*eF* 



By setting 



and 



we can write 



w*(y) := (yi - y q ci y g _i - y q c q -i , -1) , 



(5) 



tp(y,V*) = w{y*) y - y* q = w*(y) y*+y q , 

which is useful for the geometric interpretation of duality. 
Weak duality reads as follows. 



(6) 



Theorem 3.5 (|19} 125]). The following implication is true: 

[y£VAy*€V*} <p(y,y*)>0. 
Note that weak duality implies the inclusions 

V* C {y* £R q \Vy£V: tp(y, y*) > 0} and PC{i,e R q \ Vy* G V* : <p(y, y*) > 0} , 
whereas the following strong duality theorem yields even equality. 

Theorem 3.6 ([El [25]). Let the feasible sets S of ^ andT of ([D^ be non-empty. Then 



[Vy* GP* :<p(y,y*) > 0] 
[Vy£V: V (y,y*) > 0] 



ye P, 

y* G P*. 



Remark 3.7. Theorems 3.5 and \3.6[ are formulated and proved in the mentioned references 

pi 



only for the special case C 



an 



d c = (1, . . . ,1) T . However, using the generalized variants 



of scalarized problems (P\{w) ), (D\{w) ), (Pa(y) ) and (D^yjl as stated below, the general 
case can be obtained in the same way. 

The following geometric duality theorem provides a third type of duality relation. It takes 
into account the facial structure of the sets V and V* . Note that geometric duality does not 
play any role in scalar optimization because the structure of polyhedra in the objective space 
R is very simple. 
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Theorem 3.8 ([19[[TB]). VP is an inclusion reversing one-to-one map between the set of all 
K -maximal proper faces ofD* and the set of all proper faces ofV. The inverse map is given 
by 

= p| H*(y)DV*. 

ydF 

Moreover, if F* is a K -maximal proper face ofT>*, then 

dimF* +dim^(F*) = q - 1. 

Remark 3.9. The proof of the special case c = (1, 1) T and C = M? + can be found in )19$ . 
The general case can be proved in the same way using the generalized versions of ( jPi(u>)[ ), 
(Di(u;)), (P2(y)) an d (^2(2/)) as defined below. Theorem 3.8 (in the general setting) is also a 
special case of geometric duality theorem for convex vector optimization problems, see Example 
3 in JTZy. 

Remark 3.10. Non-K -maximal proper facets (faces of dimension q — 1) ofT>* correspond to 
the extreme directions ofV by a similar duality relation, where the coupling function ip has to 
be replaced by <p : M9 x R« -»• K, <p(y,y*) := (f{y,y*) + y*. The case C = R q + ,c= (1,... ,1) T 
has been studied in 125)/ and the general case is obtained likewise using the generalized variants 
of (Pi(w)), (T>i(w)), ( |P2 (j/) ) and (D^^y)) as defined below. 

The following scalarization techniques are fundamental for the algorithms described in 
the next section. As mentioned in Remarks 3.7, 3.9 and 3.10 they can also be used to prove 
weak, strong and geometric duality. The weighted sum scalarized problem for a parameter 
vector w G C + satisfying c T w = 1 is 

rp 

min w Px 



subject to Bx > b. 



Its dual problem is 



max b T u 



subject to 



B T u 



u > 



P T w 
0. 



(PiH) 



(DiH) 



The translative scalarization (or scalarization by reference variable) is based on problem 



min z subject to 



Bx > b 
Z T Px < Z T y + z ■ Z T c, 



CPa(v)) 



where Z is the matrix whose columns are the generating directions of the dual cone C + of 
the ordering cone C. The dual program is 



B T u 



max b u 



T ry 

y Zv 



subject to 

This problem can be equivalently expressed as 



P T Zv 
c T Zv 
(u,v) 



> 



(D 2 (y)) 



max b u 



T 
y w 



subject to 



B T u - P T w 

T 

c w 
Y T w 
u 



> 
> 




1 

0, 



(Da(y)) 



where Y is the matrix of generating directions of the ordering cone C . The equivalence of 
(fWy) ) and (T>2{y) ) is a consequence of the following assertion. For vectors 10 6 I', we have 



Y T w > 



Vj/GC: y T w > 



w £ C 



3v > : w = Zv. 
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4 Benson's algorithm and its dual variant 



Benson [31 H] motivated his outer approximation algorithm by practical problems that typi- 
cally involve a huge number of variables and constraints and just a few objective functions. 
He proposed three advantages of "outcome set-based approaches" in comparison to "decision 
set-based approaches". First, he observed that the set of minimal elements (in the outcome 
space R q ) has a simpler structure than the set of minimizers (in the decision space R n ) be- 
cause, typically, q <C n. This is beneficial for computational reasons but also for the decision 
maker. Second, in practice, decision makers prefer to base their decisions (at least in a first 
stage) on objectives rather than directly on a set of efficient decisions. Third, it is generic 
that many feasible points are mapped on a single image point which may lead to redundant 
calculations of " little or no use to the decision maker" [1] . 

Comparing this motivation with the notions of the previous section, we see that the solu- 
tion concepts are based on exactly the same motivation (but additionally there is a theoretical 
motivation, see |20U25| ). It is therefore not surprising that the variants of Benson's algorithm 
presented here just compute solutions in the sense of the previous section. 

The dual variant of the algorithm (based on geometric duality) has been established in 
18] . It was followed by approximating variants |31[ [32] and by a generalization of the primal 
algorithm to convex problems [9j. Unbounded problems have been first treated in |25| . 

Definition 4.1. Problem ([P]) is said to be bounded, if 

3y G R q : P[S] C {y} + C. 

The generalization from IR+ to arbitrary pointed solid polyhedral convex cones is new in 
this article but has already been used in |26j . We will present simplified variants where only 
one linear program (rather than two or three) has to be solved in one iteration^] 

The idea of the primal algorithm is to evaluate the upper image V = P[S] + C of problem 
( [P] ) by constructing appropriate cutting planes. This leads to an iterative refinement of an 
outer approximation T 5 V by a decreasing sequence of polyhedral supersets 

T° 5 T 1 5 . . . 5 T k = V. 

Both an H-representation and a V-representation of the approximating supersets T % are 
stored. The algorithm terminates after finitely many steps (say k steps) when the outer 
approximation coincides with V . 

Unbounded problems are treated by a two-phase method. First, one solves the homoge- 



neous problem (P h ) (which is unbounded, too) and its dual problem 

maximize D* h : R m x R q -)• R q with respect to < K over T (B* h ) 
with objective function 

D* : R m x R q — » R q , D*(u,w) := (wi, w q -i, 0) T . 



1 This simplification was initiated by an idea of Kevin Webster. During a lecture in the Ph.D. course in 
spring 20f f at ORFE, Princeton University, where the classical variant of Benson's algorithm was introduced, 
he proposed a variant with the two LPs ( |P2 and ( |D2(y)[ ). The advantage over the classical version is that 



( P2 (y) I and ( D2 (y) I are dual to each other. All further improvements of this article are based on this idea. 
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To this end, (P h ) is transformed into an equivalent bounded problem 



minimize P : R n -)• R q with respect to < c subject to Bx > 0, rfPx < 1, (P 7 ?) 



where 77 G int (£>* /l + K) with cr 77 = 1 (D* h denotes the lower image of (D* ft )). 



In the second phase, a primal and dual solution of the homogeneous problem (P ) are 
used to calculate a primal and dual solution of the original (inhomogeneous and unbounded) 
problem (|P]). The two-phase method requires an algorithm that works whenever an H- 
representation of an initial outer approximation 7~° 15 V with 7^, = Voo is known. If an 
H-representation of Voo is known, that is 

V 00 = {yeW\ (w i ) T y>0,i = l,...,r}, 

and if 7$ denotes the optimal value of (Pi(w/)) for i = 1, . . . , r, then 

T° = {y£Ri\ (w l ) T y> li ,i = l,...,r} 

is the desired outer approximation of V satisfying 7^ = Voo- 

If problem ([P]) is bounded, we have C = Voo-, i-e., an H-representation of the ordering C 
is required. Otherwise, whenever ([P]) is feasible, the upper image V h := P^' 1 ] + C of the 



homogeneous problem (P h ) coincides with Voo- By geometric duality, a dual solution to (P h 
yields an H-representation of V h = Voo- 

The idea of such an algorithm can be explained geometrically. Let 7~° 5 V be an initial 
outer representation of V, i.e., 7^ = Voo- First, the vertices of 7~° are computed from its 
H-representation. This can be realized by vertex enumeration, which is a standard method 
in Computational Geometry, see e.g. [2 [5]. Secondly, for a vertex t° G 7~°, problem (P2(t )) 
is solved. Usually, LP solvers yield simultaneously a solution of both the primal and the dual 
problem. If the optimal value of (P2(t )) is zero, then t° belongs to V and one proceeds with 
the next vertex of T°. If every vertex of 7~° belongs to V, we have T"° = V. Otherwise, for 



t° $ V, a solution of (P2(i )) yields a point s° G bdVPiintT , see Proposition 4.2 below. The 
solution of the dual problem, defines a supporting hyperplane H° of V that contains s°. The 
corresponding halfspace contains V but not t°. An H-representation of an improved outer 
approximation T 1 := T^nH^ is obtained immediately. This procedure is repeated until, after 
finitely many steps, T k = V. A solution (S, S h ) of ([P]) is obtained by collecting those points 
x that arise during the procedure from a solution (x,z) of (P2(i 1 )) with zero optimal value. 
In this case we have t = Px for some vertex t of T l . Hence Px is a vertex of V which implies 
that x is a minimizer for (|P]). In the unbounded case, S h contains directions that originate 
from a solution to the homogeneous problem. A solution of the dual vector optimization 



problem (D* ) is obtained by collecting those dual solutions of (P2(t J )) with non-zero optimal 
value. 

Proposition 4.2. Let S 7^ 0, C C M q a solid pointed polyhedral cone and let c G int C . 

For every t G U. q , there exist optimal solutions (x,z) to (P2(i)) and (u,w) to (D2(i)). Each 
solution (u,w) to (D2(i)) defines a supporting hyperplane H := {y G R q \ w T y = b T u] ofV 
such that s := t + z ■ c £ H DV . We have 

t^V z>0, tGwMinP z = 0, t £ intV z < 0. 
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Proof. Fix t G M 9 . Since 5/8 and c G intC, (P 2 (i)) is feasible. Assuming (P 2 (*)) is 
not bounded, we obtain t + (z — n)c — Px G C for all n G N. Dividing by n and letting 
n — > oo, we conclude — c G C. As c G int C, convexity of C implies G int C. Thus C = M 9 , 
a contradiction. Consequently, (P 2 (t)) and, by duality, also (D2(t)) have optimal solutions 
(x, z) and (u, u)), respectively. Duality yields b T u — t T w = z and thus s = t + zc belongs to 
H. Of course, H is a hyperplane as the constraint w T c = 1 of (D2(i)) implies u) 7^ 0. From 
Px <c t + z ■ c we conclude that s = i + zc belongs to V. For arbitrary y G "P, there exists 
x G 5 such that 1/ >c Px. Hence (x,0) is feasible for (P2(y)). Weak duality between (P2(y)) 
and (D2(y)) implies that y T w > 6 T u for every (u,w) G T, in particular, t/ t, u) > Hence 
H = {y G M 9 ! = b T u} is a supporting hyperplane to V . The remaining statements are 
now obvious, where the fact wMin"P = bd"P can be used. □ 

Proposition 4.3. Every vertex ofV is minimal. 

Proof. Let y G V = P[S]+C be not minimal for V . Then there are z£P and k G C\{0} such 
that z = y — k. The points y — k and y + k belong to V and we have y = ^(y — k) + ^(y + k). 
Hence y is not a vertex of "P. □ 

Two functions are used in the following algorithm. The function dual() computes a V- 
representation of an outer approximation T from an H-representation of T, i.e., this function 
consists essentially of vertex enumeration. This H-representation of 7~, however, is stored in a 
dual format, namely, as a V- representation of an inner approximation T* of the lower image 



V* of (|pj), where 

T* = {y* eM. q \<p(y,y*)>0,yeT}. (7) 
The following duality relation holds. 

Proposition 4.4. If 7^ T C M q is closed and convex and Too 5 C, then 

T={yGM' ? | <p(y,y*) >0,y* G T*} • (8) 

Proof. The inclusion C is obvious. Assume that the inclusion D does not hold, i.e., there 
exists y £ M q \T with (p(y,y*) > for all y* G T* ■ Applying separation arguments, we 
get r\ G C + \ {0} with rj T y < mi y& -j-r] T y := 7. By ([3]), we can assume rj T c = 1. Setting 
y* := (771,.. . ,t? 9 _i,7), we get = rj and <p(y,y*) = w(y*) - y* = ify-^y. For all y G 7", 

we have r/ T y — 7 > 0, i.e., y* G 7"*. But </?(y, y*) = r/ T y — 7 < 0, a contradiction. □ 

In the following algorithm, a V-representation of a polyhedron T (that contains no lines) 
is denoted by (T po ints,Tdirs), i-e., T = convT po ints + coneTdirs- We assume further that a V- 
representation returned by the function dual() is minimal, i.e., Tpoints consists of only vertices 
of T and 7dirs consists of only extreme directions of T. 

The function solve() returns an optimal solution (x, z) of (P2(£)) and an optimal solution 
(u, w) of (D2(t)). Since (D2(i)) is, up to a substitution, the dual program of (P 2 (i)), only one 
linear program has to be solved. 

The variables in the following algorithm are arrays of vectors. By |^4| we denote the 
number of vectors in an array A and by A[i] we refer to the i-th vector in A. The command 
break exits the inner-most loop. 
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Algorithm 1. 



Input: 

B,b,P,Z (data of Q); 

a solution ({0} , S h ) to jF}; 

a solution f h to flD* fc | ); 
Output: 

... a solution to Q; 
f ... a solution to ( |D*j ); 
(T p0 ints,Tdirs) • •• a V-representation of P; 
(T^i n t s , {~e q }) ••• a V-representation ofV*; 

f <- {(solve(Di(u;)),u/)| («,«;) G f h }; 
repeat 

/Zag <- 0; 

5^0; 

(%oints,T di rs) <- dual(7^ nts ,{-e«}); 
for i = 1 to |7^omts| do 

^ ^ Tpoints [ ^ ] > 

(x,z,u,«;) <- solve(P 2 (t)/D 2 (i)); 
if z > then 

^TU{( UlW )}; 

/Za# <- 1; 

break; (optional) 
else 

S ^ Su{x}; 
end if; 
end for; 
until flag = 0; 



Theorem 4.5. Ze£ S suppose P h has a vertex and assume that the command 

{Tpoints-, Tdirs) ^~ dual (7^ j„j s , Tdirs) 

generates a minimal V-representation of T from a given V-representation of T* according to 
([7]). Then Algorithm 1 is correct and finite. 

Proof. As T h is non-empty (by the definition of a solution), we can choose some (u, w) 6 T h . 
Then D* h (u,w) = (w\, . . . ,w q —i,0) is if -maximal in T>* h . Hence u solves the homogeneous 



variant (i.e., we set b = 0) of {Di(w) ). Consequently, (T>i(w)) (for arbitrary b) is feasible. 



Since 5/0, (Pi (it;)) is feasible, too. Thus, by linear programming duality, (Di(w)) has a 
solution. 

The set T* := conv7I,gj n ^ s + cone { — e^} is a non-empty subset of T)* . Hence, by Theorem 



3.5, after calling the function dual(), T := conv T po ints + coneTdirs is a superset of V . 

As T h solves the dual of the homogeneous problem, we have Too = Poo = P h , see 
Section 4.6] for more details. As V is assumed to have a vertex, T must have a vertex, hence 
the array Tpoints is non-empty. 
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By Proposition 4.2 solutions to (P2(i)) and (T>2(t)) exist. The vectors x £ S are mini- 
mizers of Q. Indeed, x is added to S only if z = 0. In this case, we have t EV, where t is 
a vertex of T 5 because, by assumption, T po ints contains only vertices of T ■ Hence t is a 
vertex of P and, by Proposition 4.3 i is a minimizer for (|P|). 

The algorithm terminates if all vertices of T belong to V . Since Voa = 7^o, we conclude 
V = T, i.e., (S, S h ) is an infimizer of ([P]) and {T po ints,Tdirs) is a V-representation of "P. 

i.e., T consists of only maximiz- 
es. "' ' " " 



A solution (u, w) to (D2(t)) is always a maximizer of (D* 
Since at termination T 
for ©and (T p * omts ,{ 



V, Theorem 3.6 implies T* 
-e q }) is a V-representation of T>* . 



T>* and thus T is a supremizer 



Finally we show that the algorithm terminates after a finite number of steps. The point 
:= t k + z k ■ c computed in iteration k (consider the 'repeat' loop) by solving (P2(t fc )/D2(i fc )) 



belongs to int T k 
by Proposition 



4.2 



1 whenever z k > 0. 
we know that F :- 



We have T k : 
- {y G V\ (w^ 



r fe_1 n {y G M 9 | (w k ) T y > (u fc ) T 6} and 
= (u k ) T b} is a face of V with s fc G F, 



where F C bd7 _/c . This means for the next iteration that s k+1 i 7 (because s k+l G intT" fc 
and therefore s k+l belongs to another face of V . Since V is polyhedral, it has a finite number 
of faces, hence the algorithm is finite. □ 

We now turn to the dual variant of Algorithm 1 . An analogous construction is now applied 
to the lower image T>* , i.e., a finite sequence of polyhedral sets 



*1 



D T* 



V* 



is calculated. Using the upper image V h (which is a polyhedral cone) of the homogeneous 
problem (PN, we define the set 



A := {y* G R q \ w(y*) G (V h ) + } 



The counterpart of Proposition 4.2 is the following. 



Proposition 4.6. Let S ^ and t* G A. For w := w(t*), (Pi(w)) has a solution and for 
every such solution x, H* (Px) is a supporting hyperplane of T>* that contains 



(t* v ...,t* 1 ,w T Px) G Max K V*. 



(9) 



Moreover, we have 
t* G" V* 



w 



T Px < t* 



t* G Max^P* 



w Px 



Proof Since t* G A, for all k G V h , we have w T k > 0. This means that the homogeneous 
variant of the linear program ( |Pi(u>)[ ) (i.e., we set b = in (Pi(w) )) is bounded (and feasible, 
as is feasible). Consequently, the dual program is feasible, even for arbitrary b, i.e., ( Di(u>)| ) 
is feasible. On the other hand, (Pi(w)) is feasible, since we assumed S ^ 0. Altogether this 



implies that both (Pi(w)) and (Di(w)) have optimal solutions denoted, respectively, by x 
and u. Strong duality implies w 1 Px = b T u. Thus, ([9J holds. We have s* G H*{Px) because 
this can be written as w(s*) T Px = s* where we have w = w(t*) = w(s*). Together with 



Theorem 3.5 we obtain that H*{Px) is a supporting hyperplane of T>* that contains s* . The 
remaining statements are now obvious. □ 



The following consequence of Proposition 4.4 is useful to characterize the condition t* G A. 
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Corollary 4.7. Let the assumptions of Proposition 4-4 be satisfied. Then, w{y*) G (7^o) + 
for all y* eT*. 

Proof. Assuming the contrary, there is y* G T* and k G Xio with w(y*) T k < 0. Let y G 7". For 
sufficiently large A > 0, using ([6]), we obtain tp(y + \k,y*) < 0, which contradicts Proposition 

□ 



The following dual algorithm has the same input and output as Algorithm 1 . Similar func- 
tions are used. The function dual() computes a V-representation of an outer approximation 
7~* of T>* from a V-representation of an inner approximation 7~ of V . In contrast to Algorithm 
1, it is not necessary that dualQ returns a minimal V-representation. The recession cone of 
sets 7~* occurring in the algorithm is known, in fact, we always have 7^, = —K = M.+ {— e q }. 
Therefore we denote the return of the function dualQ by (T* oints , ~) indicating that the second 
return value (the array containing the extreme directions of 7~*) is not used. 

The function solve () returns an optimal solution x of (Pi (to)) and an optimal solution u 
of (Di(u;)). Again, only one linear program has to be solved. 



Algorithm 2. 



Input: 

B,b,P,Y (data of Problem Q); 
a solution ({0} , S h ) to (Q; 



a solution T h to flD* 7 ' 
Output: 

(S,S h ) is a solution to Q; 



T is a solution to (D* 
(T p0 ints,Tdirs) ••• a V-representation of V; 
0~p* oints , i~ e q}) ••• a V-representation of V*\ 



{Px\ x G S h }u{y\ y is a column of y}; 
{w| (u,u>) G f h }; 



Tdirs 

w <— mean 
S <- {solve(P 1 (u)))}; 
repeat 

flag <- 0; 

Tpoints 4 {-Px| X G S*} 
(joints' ~) ^ dual (^points, Td 



1 to 



'points 



Tpoints ' 



do 



for i 

(a;, it) <(— solve(Pi (w)/Di (tu)); 



if t* - b T u > then 

5^ Su{z}; 
/Zag <- 1; 
break; (optional) 
else 

T<-T U{(u,w)}; 
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end if; 
end for; 
until flag = 0; 

delete points x E S whenever Px is not a vertex of T ; 



Remark 4.8. The last line in the algorithm is easy to realize, for instance, by computing a 
minimal V -representation using the command 

{Tpoints j Tdirs) dual {Tp i n tsi Tii rs ) 

from Algorithm 1 by standard vertex enumeration methods. Then one has to test if for x E S, 
Px belongs to T po ints, if not, x is deleted from S. In particular, it is not necessary to solve a 
linear program. 

Theorem 4.9. Let S ^ and assume that V h has a vertex. Then, Algorithm 2 is correct 
and finite. 



Proof. By similar arguments as in the proof of Theorem 4.5 one can show that Pi(iu)) has a 
solution. 



The set T := conv Tpoints cone Tiirs is a subset of P. Hence, by Theorem |3 . 5| , after 
calling the function dual(), T* := convT* D j n4s + coneTJ^,, is a superset of V* . Since T ^ M 9 , 
c E intC and To =! C, we have cone 7^ = M+ • {— e q }, i.e., we can set T^ irs = {— e q } and we 
know that T^ nta ^ 0. 

The array ({0} ,Tdi rs ) provides a V-representation of P h , i.e., To = P h ■ Corollary 



yields that T* oints Q A. Hence, by Proposition 4.6 solutions to (Pi(w;)) and (Di(io)) exist. It 



4.7 



can be easily shown that the vectors (u,w) E T are maximizers of (D* ), see also [25l Lemma 
4.51]. 

The algorithm terminates, if T* oints C V* . Since V*^ = T£ = R + {-e q }, we conclude 
V* = T* , i.e., T is a supremizer of (D* ) and {T* oints , {— e q }) is a V-representation of V* . Since 



at termination T* =V*, Theorem 3.6 implies T = P. Thus {Tpoints-, Ti r ) is a (not necessarily 
minimal) V-representation of P. A solution x to (Pi(u;)) is in general not a minimizer for 
( |P| ) (but only "weakly efficient", compare e.g. |25l Theorem 4.1]). Therefore, in the last 
line of the algorithm, x is deleted from S, whenever Px is not a vertex of T ■ According to 



Proposition 4.3, the remaining set S consists of only minimizers. It is non-empty because, 
by assumption, P h has a vertex and hence T = P must have a vertex. As non-vertex points 
are redundant in a V-representation of a set which has a vertex, the property of S being an 
infimizer for (|P]) is maintained by deleting the non-minimizers in S. 

Finally we show that the algorithm terminates after finitely many steps. We consider 
the 'repeat' loop in iteration k. We set w = w k and t* = t* k and denote the solutions 
of (Pi(w k )) and (Bi(w k )) by x k and u k , respectively. The point s* k := t* k + z* k ■ {-e q }, 
where z* k := (t* k — b T u k ) belongs to p*i k ~ l ) \ Max^ p*^- 1 ) whenever z* k > 0. We have 



4.6 



we know that F* 



T* k := T* (fc_1) n {y* E R q \ Lp{Px k ,y*) > 0} and by Proposition 
{y* E T>* | ip(Px k ,y*) = 0} is a face of T>* with s* k E F*. Likewise to [25, Lemma 4.48], 
we see that F* C MaxxT* k . This means for the next iteration that s*( fc+1 ) ^ F* (because 
s *(fe+l) e T *k ^ Max K T* k ), and therefore s *( k + l ) belongs to another face of V*. Since V* is 
polyhedral, it has a finite number of faces, hence the algorithm is finite. □ 
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Let us summarize the two-phase method for solving unbounded problems. We consider an 
arbitrary linear vector optimization problem, where we only assume that C is a solid pointed 
polyhedral cone. We fix some c according to ([3]), which is always possible in the way described 
after (|3j. In phase 1, we first try to compute some rj G int (T>* + K) with n T c = 1. This can 
be realized by Algorithm 3 in |25} Section 5.5], where the set T has to be adapted to the more 
general setting of this article. Note that c has a different meaning in \25\ Section 5.5]. The 
first LP solved by the mentioned algorithm is 



rp rp 

min w + u 



s.t. (u, w) G T. 



If this linear program is infeasible, then (D* ) is infeasible. Otherwise one obtains either some 



rj G int (T>* + K) or the algorithm indicates that int (T>* + K) is empty. In the latter case, we 
know that V h has no vertex. This means that V, if non-empty, contains a line. This case has 
not been treated so far. Since c q = 1 according to ([3]), the condition n T c = 1 can be always 
realized by an appropriate choice of r/ q . 



Next, we solve (P'') by either Algorithm 1 or Algorithm 2. Since (P^) is bounded, a 



solution of the primal and dual homogeneous problem of (P v ) can be easily obtained. However, 



this not necessary as the 14-components of (u, w) G T h are not used in Algorithms 1 and 2. 
Therefore we can use 



S" 



and 



0. 



z is a column of Z 



as an input of Algorithm 1 or 2 to solve (P v ), compare also [251 Theorem 5.20]. Let (S^, 5") 
be a solution of and let T v be a solution of the dual problem of Then, a solution 



of {P h ) is obtained by setting 



(5,5") := ({0},S„\{0}), 



compare [23 Theorem 5.23]. A solution T h of §D*H can be obtained from 7^ but again only 
the ui-components are required by Algorithms 1 or 2 in phase 2. As a consequence of [251 
Theorem 5.25], we can use 

f h :={(0,w(y*))\y*eT; oints , y* q = 0} 

as an input of Algorithm 1 or 2 in the second phase, where 7~* oints is the result from the 
first phase, i.e., (7^, ints , {— e 9 }) is a V-representation of the lower image V* v of the dual 



problem of (P 77 ). In the second phase, the first LP to be solved in Algorithm 1 is (T)i(w)). 



If (Di(u>)) turns out to be unbounded, we know that (|P]) is infeasible. Likewise, if the first 
LP in Algorithm 2, namely (Pi(tt))), is infeasible, we know that (|P]) is infeasible. Otherwise, 
according to Theorems 4.5 and |4.9[ solutions of ([P]) and (D* ) are computed. 



Remark 4.10. In practice the condition z > in Algorithm 1 is replaced by z > e for some 
e > 0. Consequently, Algorithm 1 yields an e-solution of ([P]) in the sense that in Definition 
3.2 the finite infimizer is replaced by a finite e-infimizer, i.e., condition ^ is replaced by 



conv P[S] + cone P[S h ] + C - e {c} D P[S] + C. 



(10) 



Taking into account that (using the assumption c q = 1 in (|3 

<f(y - sc, y*) = ip(y, y* + ee q ), 
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we see that Algorithm 1 yields an e-solution of (D*) in the sense that in Definition 3.4 



finite supremizer is replaced by a finite e-supremizer, i.e., condition Q is replaced by 

conv£>*[f] -K + e{e q } D D*[T] - K. (11) 

Likewise, in Algorithm 2 the condition t* — b T u > is replaced by t* — b T u > e for some 
e > 0. Consequently, Algorithm 2 yields an e-solution of ([P]) and an e-solution of (D* ). Note 
that an e-solution (S, S h ) of ([P]) refers to an inner and an outer approximation of the upper 
image V in the sense that 

conv P[S] + cone P[S h ] + CQVQ conv P[S] + cone P[S h ] + C - e {c} . 



Likewise, an e-solution T of (D* ) refers to an inner and an outer approximation of the lower 
image T>* in the sense that 

conv D* [f] - K C V* C conv D*[f] — K + e {e q } . 

Note that the approximation error of the classical variant of Benson's algorithm and its dual 
variant has been studied in 1311 [31 



5 Computation of polyhedral set-valued risk measures 

Set- valued risk measures evaluate the risk of multi-variate random portfolios X : — > ~K d the 
components Xi (w) of which represent the number of units of the i-th asset in the portfolio, 
i = 1, . . . , d. If transaction costs are present, such risk measures are more appropriate than 
real- valued functions, which always represent a complete risk preference and thus cannot 
account for incomparable portfolios. 

The theory of set-valued risk measures was initiated in [23] and systematically developed 
in |12j and [H] . We refer the reader to these references for further motivation and information. 
Here, we restrict ourself to the case of finite probability spaces and the question how the values 
of a set-valued risk measure can be computed. It will turn out that this leads to problems of 
type (P), hence one can apply the algorithm presented in Section |4j 

The basic idea is as follows. The value of a set-valued risk measure at some random future 
portfolio X consists of initial deterministic portfolios u which can be given as deposits for 
the 'risky payoff' X , thus making the overall position 'risky payoff plus deposit' a non-risky 
one. It usually is not possible to use all assets as deposits, but rather a small subset including 
cash in a few currencies, bonds, gold or similar risk-free or low-risk assets. These 'eligible' 
assets are assumed to span the linear subspace M of M. d with 1 < dim M = m < d. A typical 
example, already used in [23J, is W n x {0} d m , i.e. the first m assets are eligible. 

Let (fi, P) be a finite probability space and N > 2 be the number of elements in Q = 
{wi, . . . , wjv}- We assume p n = P ({u n }) > for all n £ {1, ...,JV}. The space of all 
multi-variate random variables X : f2 — > M. d is denoted by L^(Q,P). A random variable 
X 6 L? d (Q, P) can be identified with an element x £ ]& dN through 

x = (X 1 (u 1 ),...,X d (io 1 ),X 1 (oo 2 ),...,X d (io N )) T eR dN 

and vice versa. Thus, the function T: L d (Q,P) — > M. dN defined by TX = x is a linear 
bijection. If A C L° (O, P) then we set A = [x G R dN \ x = TX, lei}. 



16 



Let Kq C M be a finitely generated convex cone satisfying W\_ C Kq ^ Mr. Such a 
'solvency' cone models the market at initial time. We set Kq 1 = Kq n M and P (M, 1£q ) = 
{DCM|D = D + K$ { }. A risk measure is a function P: L° (ft, P) V (M, K$f) satisfy- 
ing 

Vu G M, VX G (ft, P) : P (X + ul) = R (X) - u (12) 
R with 1 (w) = 1 for all w G ft is the uni-variate random variable with constant 



where 1: ft — 
value 1. 

With P, we associate a risk measure P: 
for x = TX. Consequently, P satisfies 



V (M, Kq 1 ) by means of R(x) = R (X) 



Vu G M, Vx G M d7V : P ( x + I d u ) = P (x) 



1i 



(13) 



where 



G 



ndNxd 



and Irf 



/ 1 






1 



... \ 
... 

1/ 



ndxd 



The most common way to generate a risk measure is by means of a set A C (ft, P) of 
random variables which are considered to be 'acceptable' by the decision maker. The value of 
a risk measure generated by A then consists of all deterministic (available at time 0) portfolios 
u G M C Mr which, when added to the uncertain future position X, make the overall position 
acceptable. Thus, 

R A {X) = {u G M | X + ul G A} . 



This functions indeed satisfies (12). Correspondingly, 

R A ( x ) = ju G M | x + I d u G A\ 

satisfies ( fl3| ), and we have R^ = Ra- Vice versa, with risk measures P: L°(ft,P) — > 
V (M, K^) and R: R dN ->• P (M, P^ f ) we associate the sets 

A R ={X £ L° d (Q,P) | G R(X)} andA R = jx G R dN | G P(x)} , 

respectively. A basic fact about risk measures is a one-to-one correspondence between closed 
acceptance sets A Q L d (ft, P) which satisfy A + K^ 1 ! C A and risk measures P: (ft, P) — )• 
P (M, ifg ) with a closed graph by means of the above formulas. In particular, the relation- 
ships A = Ar. and P = Ra^ hold true. See [Hj for further details. 

A risk measure P is called polyhedral if the associated risk measure P is polyhedral, i.e., 

if 



gr P 



u) G 



t>dN 



x M I u G P (x 



»} 



is a polyhedral subset of M dN x M d . The one-to-one correspondence between risk measures 
and their acceptance sets extends to the polyhedral case: A C M dN is polyhedral if, and only 
if, P^ is polyhedral, and P is polyhedral if and only if Ar is polyhedral. 
The above discussion leads to the following conclusion. 
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Remark 5.1. Since each polyhedral risk measure R has the representation 

R(x) = |ii g M | x + I d u g Afi\ 

where Ap, is a polyhedral set, the set R (x) is the upper image of a linear vector optimization 
problem. Indeed, if Ap, C M dN has the H-representation Ap, = |y G M dN \ By > &| where B 
and b are matrices of appropriate dimension then 

R(x) = ju G M | B (x + i d uj > S} = |u G M \ BI d u >b- Bx} . 

Let P G ^ dxm 5 e a matrix with column vectors /z 1 , . . . ,/x m forming a basis of M and define 
B = BI d P, b = b — Bx. Then, observing that R(x) + Kq = R(x) and substituting u = Pz 
we obtain that R (x) is the upper image of the problem 

minimize P : M m — > M with respect to ^k^ 1 SU0 3 ec t to Bz > b. 

However, this is just a theoretical result since in practice life is not as straightforward: Usually, 
the constraints describing A involve a large number of auxiliary variables, and u is given as 



a linear function of those (see Example 5.3 below). Therefore, the algorithm presented in 
Section^ is an appropriate tool to compute the values of a polyhedral set-valued risk measure 
because the dimension of the pre-image space usually is much greater then the dimension of 



the image space which is dimM = m < d. Compare Examples 5.2 and 5.3 below. 



In the following, we will discuss two examples which will be used for the numerical com- 
putations reported in Section [6| 

Example 5.2. In worst case, the regulator/decision maker only accepts positions with 
non- negative components. Thus, the acceptance set is A = (-^)+ which is the set of all 
component- wise non- negative random variables. The market extension of the worst case risk 
measure, i.e. when trading is allowed, is related to the set of superhedging portfolios, see [26] . 
Its acceptance set in a one-period market is A = (L d ) + + A'ol + L d {Kt) where the cone Kq 
and the random cone Kj> model market conditions with a potential bid ask price spread at 
initial and terminal time, respectively, and L° d (Kt) = {X G L d \ G Q : X (w) G Kt 
The cones Kt (w) are also finitely generated convex cones satisfying C Kt (w) ^ M rf for 
all co G fi. 

The market extension of the worst case risk measure still is very conservative since a 
payoff X is acceptable only if there is a trading strategy such that its result, when added to 
X, is non- negative in all components in all possible scenarios, even those with a very small 
probability. Therefore, we introduce a "relaxed" variant as follows. 

Let G C M. d be a finitely generated convex cone with C G ^ R d and consider the 
following acceptance set 

a rwc = | x g L° d (n,P) | Vw G n-. X(uj) G (-e + K^) n G} + K 1 + L° d (K T ) 

where e G M d such that > for all i G {1, . . . , d}. Compared to the " true" worst case risk 
measure, the set (L d ) + is replaced by 'a little' bigger set. 

Thus, payoffs with "small" negative components may still be considered acceptable, and 
the size of the risk related to such payoffs is controlled by e and G. The cone G may serve 
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as a conservative estimate of a market model which the regulator /supervisor thinks is robust 
enough to cover most market scenarios. For example, G can be chosen such that P{G C 
Kt) > 1 — a for some significance level a G [0, 1]. Then, the probability of a loss is bounded 
by a, and a potential loss (in physical units) is bounded by e. Note that in the scalar case 
d = m = 1, the relaxed worst case risk measure reduces to the scalar worst case risk measure. 
The market extension of the relaxed worst case risk measure is given by 

RWC (X) ={u GM\zeK , Z eL° d (K T ) ,VwGO: 

i(w)-z-Z(w) + «e(- f + if)nG} (14) 

and can be seen as a relaxation of the superhedging set and thus as a good deal price bound 



of —X, see Example 6.6 for details. It is polyhedral (convex) as A RWC is polyhedral (convex), 
but not sublinear. This is a new feature since the classical worst case risk measure is always 
sublinear. 

In order to describe RWC (x), let g l , . . . , g be the generating vectors of the cone G and 
let G be the dN x LN matrix which contains N blocks on its diagonal, where each block 
consists of the matrix with g 1 , . . . , g L as columns. Then 

Vn € {1, . . . , N} : X {uj n ) G G ^ 3 7 G R+ N : x = &y. 

Similarly, let Kt be the dN x J matrix which contains N blocks on its diagonal, where 
the first block consists of the matrix with the generating vectors of Kt{<jJi) as columns, the 
second block contains the generating vectors of Kt(co2) and so forth to the last block with 
the generating vectors of Kt(ojn)- J is the sum of the number of generating vectors of all 
Kts. Let Kq denote the matrix containing the I generating vectors of Kq as columns. Then, 

RWC (x) ={Pz\x + I d (Pz + e- K r^j - K T s G R+ N , 

i d (Pz - K r^j - K T s -G"/ = -x,z£ R m , 7 G R+ N , r G s G R{}, 

with ordering cone Kq 1 . Thus, the dimension of the pre-image space is m + LN + / + J 
whereas the dimension of the image space is just m. 

Example 5.3. The following set-valued function is a generalization of the scalar average value 
at risk (see pUl p. 210]) which is probably the most important and most studied example of 
a sublinear coherent measure of risk as introduced by pQ. Let a G (0, l] d . Define for X G L d 

AV@R a {X) ={diag (a)" 1 E [Z] - z \ Z G (L° d ) + , 

X + Z - zl G K 1 + L° d (K T ) , z G R d } Pi M (15) 

where diag (a) -1 is the inverse of the diagonal matrix with the components of a on its main 
diagonal and zero elsewhere, and the cones Kq, Kt (w) modeling the market conditions are as 
described above. Therefore this risk measure is also called the market extension of a simpler 
'regulator' version, see |15j . We also refer to this paper for further motivation, interpretation 
and more details. It is immediately clear that AV@R a is not given in the form of Ra above, 
but it is a polyhedral convex (even sublinear) risk measure. 
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Its 'hat' variant can be derived as follows. Replace Z by z £ Mr and introduce auxiliary 
variables which admit to write x + z — I^z as non- negative linear combination of the generating 
vectors of the cones Kq and Kt (w n ). Transform the objective into matrix form and get 

AV@R a (z) = {Px \Bx>b} + K^, 

with appropriate matrices B, b and P where = Kq n M as before. The dimension of the 
pre-image space is d(N + 1) + I + J whereas the dimension of the image space is just m. It 
is worth mentioning that in the scalar case (i.e., without transaction costs) Rockafellar and 
Uryasev observed that the AV@R can be computed by solving a linear optimization problem, 
see [29]. 



6 Numerical examples 

The algorithms have been implemented with MATLAB using the GNU Linear Programming 
Kit (GLPK) to solve the LPs and the CDDLIB package [5] for vertex enumeration. The 
graphics have been generated by Java View^] and OpenOmce (Figure [H]). By a straightforward 
extension of the above results we can also solve linear vector optimization problems with 
constraints of the form 

a<Bx<b lb<x<ub, (16) 

where the components of a, b, lb, ub belong to RU {— oo, +oo}. All examples were computed 
on a MacBook Pro with 2.26 GHz clock and 8 GB memory. We made use of the fact that 
all the LPs have a very similar form. This means that the matrix B does not need to be 
changed during the algorithm (except one line because r\ is not yet known at the beginning). 
This allows us to initialize LPs by appropriate basis solution of LPs solved in previous steps 
(warm starts). 

In the following examples we provide tables with a few computational data, such as the 
total time and the number of LPs solved LPs). Note that we compute an e-solution (S, S h ) 



of (|Pj) and an e-solution T of (D* ), compare Remark 4.10 We provide the cardinality | . | of 
the sets S, S h and T. Recall that we have \S h \ =0 whenever the problem is bounded. Note 
that \S\ and \T\ "correlate" to the number of, respectively, vertices and facets of V (but do 
not need to coincide exactly). One reason for possible differences is degeneracy as discussed 
in |25t Section 5.6], another one is numerical inaccuracy. 

Further we denote by t max the maximum time used to solve one LP and by t aver the 
average time to solve one LP. The quotient t max /t aver indicates the impact of using warm 
starts. We start with two numerical examples from the literature. 



Example 6.1. Shao and Ehrgott [3T] used extended variants of Benson's algorithm to solve 
linear vector optimization problems occurring in radio therapy treatment planning. We com- 
pute Example (PL) in |3l] which has three objectives and a matrix B of size 1211 x 1143 with 
153 930 nonzeros. The ordering cone is C = M$, and the problem is known to be bounded, 
which means that the first phase of our algorithms as well as the computation of r/ can be 
skipped. Further we set c = (1, 1, 1) T . 

The following table shows some results obtained by Algorithm 1. The second column in 
the table concerns the optional break command in the algorithm. One can observe that more 



by Konrad Polthier, http://www.javaview.de 
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Figure 1: Illustration of the upper image V in Example |6.1| top left: primal algorithm, variant 
'break', e = 0.3; top right: primal algorithm, variant 'no break', e = 0.3; bottom left: dual 
algorithm, variant 'break', e = 0.3; bottom right: primal algorithm, variant 'break', e = 0.05. 

LPs have to be solved when the break command is disabled. On the other hand, less vertex 
enumerations are required. This explains why the variant 'no break' is becoming faster than 
the 'break' variant when e > is small enough. 



£ 


variant 


total time 


\s\ 


\T\ 


#LPs 


tmax 


t-max /taver 


0.3 


break 


47 sees 


46 


29 


75 


0.84 sees 


1.8 


0.1 


break 


91 sees 


104 


61 


165 


0.87 sees 


2.0 


0.05 


break 


144 sees 


176 


94 


270 


0.86 sees 


2.0 


0.005 


break 


1596 sees 


1456 


597 


2053 


0.84 sees 


1.9 


0.3 


no break 


54 sees 


54 


34 


88 


0.85 sees 


1.8 


0.1 


no break 


114 sees 


134 


78 


212 


0.84 sees 


1.9 


0.05 


no break 


205 sees 


264 


129 


393 


0.85 sees 


1.9 


0.005 


no break 


1411 sees 


1945 


804 


2749 


0.84 sees 


1.9 



Although we need less computational time than in [31], it is difficult to compare the 
results as we use a faster computer, a different (open source) LP solver, and we utilize warm 
starts. Moreover, in [31] an online vertex enumeration method is used, which is preferable 
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if the number of vertices and facets of V is large. Furthermore, our method yields the same 
approximation error as in |31| by less vertices and facets of V . See Figure [I] for an illustration 
of part of the results. 





Figure 2: Illustration of the upper image V in Example 6.2 for ej ||c|| 
approximation, left) and for ej ||c|| = 10 -6 (right). 



10 (inner and outer 



Example 6.2. Ruszczyhski and Vanderbei [30] developed a specialized parametric method for 
computing all minimizers of bi-criteria problems. Using intermediate results of the parametric 
simplex method, they solved in |30j . for instance, a mean-risk model with a dense matrix B 
of size 6161 x 3799 having 4 435 919 nonzero entries. They pointed out that computing all the 
5017 minimizers takes only a little more time than solving one single LP. As this problem is 
known to be bounded, we can skip in our algorithms the first phase as well as the computation 
of 77. For c = (1, 1) T , our primal algorithm yields approximate solutions as shown in the 
following table. 



e 


total time 


\S\ 


\T\ 


#LPs 


tmax 


tmax /taver 


y/2- 10- 4 


946 sees 


6 


7 


13 


347 sees 


4.4 


V2 ■ 10~ 5 


1648 sees 


22 


23 


45 


304 sees 


8.9 


V2 ■ 10" 6 


3085 sees 


62 


63 


125 


310 sees 


14.1 



We see that a "good" approximation with ej ||c|| = 10 6 can be obtained in about ten times 
the time required to solve a single LP. This means our general method needs much more time 



for an e-solution (compare Remark 4.10) than the parametric method for bounded bi-criteria 
problems in [30] needs for the exact solution. On the one hand, approximating solutions are 
often sufficient for a decision maker in practice, compare Figure [2} On the other hand, we 
think that the ideas of the algorithm by Ruszczyhski and Vanderbei are promising for further 
improvements of Benson type algorithms for arbitrary linear vector optimization problems. 



The following three numerical examples refer to Example 5.3 in the previous section. 



Example 6.3. Let us consider d = 12 assets, the first one is a risk-free USD bond with 
annual interest rate 5%. Given is the vector of today's asset prices, the vector of the expected 
returns and a covariance matrix for the other 11 correlated risky assets denoted in USD. 
Then, one can set up a one-period tree for the asset prices St with time horizon T = 1 year 
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0.675 0.7 0.725 0.75 0.775 



Figure 3: Illustration of the upper image V (left) and the lower image T>* (right) in Example 
Ufor e = 10~ 4 . 



as in [21] to reflect the drift and covariance structure. The resulting number of scenarios is 
N = 2 d ~ 1 = 2048. We consider proportional transaction costs for the bond to be Ao = 3% and 
for the first risky assets (usually another currency) to be Ai = 7%, the second risky asset to 
be A2 = 5% and all other risky assets to be 1%. Then, the bid and ask prices of the assets are 
(Sf)i = (1 + Xi)(S t )i and (sf)i = (1 - Xi)(S t )i for % = 0, . . . , 11 and t € {0,T}. Furthermore, 
let us assume an exchange between any two risky assets can not be made directly, only via 
cash in USD by selling one asset and buying the other. Since the risk-free bond has strictly 
positive transaction costs Ao, the cones Kq and i£r(^n) for n = 1, . . . , N have d(d — 1) = 132 
generating vectors each. Thus / = 132 and J = 270 336. 

We want to evaluate the risk of an outperformance option with physical delivery and 
maturity T. This option gives the right to buy the asset that performed best out of a basket 
of assets at a given strike price. Let the strike be K = (1 + Ai)(<Sq)i. To normalize to today's 
prices, let a vector g be defined as (Sq)i = ^(S'o)^ for i G {1, . . . , 11}. The payoff X of the 
option is — K in the risk free asset, gi units of asset i for the smallest i satisfying gi(S^) 1 = 



max 



.,ii}(9j(ST, 



3 ) > K and zero in the other assets. If max 



J6{1, 



,ii}(9j(Sj<. 



J ) < K the 



payoff is the zero vector. 

Let us calculate AV@R a (X) as described in Example 5.3 with significance levels 



a 



(0.1, 0.08, 0.09, 0.1, 0.05, 0.05, 0.04, 0.05, 0.03, 0.04, 0.03, 0.04) rj 



As the space of eligible assets we choose the space spanned by the first and the second asset, 
i.e. M = M 2 x {0} 10 . Formula (15) leads to a linear vector optimization problem with 2 
objectives and constraints of the form (16) where the matrix B is of size 24 586 x 295 056. B 
is sparse having 1 150 986 nonzero entries. The ordering cone is Kq , which is strictly larger 
than and is generated by 2 vectors. The vertices of AV@R a (X) are minimal deposits in 
the bond and the second asset that compensate for the risk of X measured by AV@R a . The 
following table shows some computational data of the primal algorithm. 



£ 


total time 


\S\ 


\S h \ 


\T\ 


#LPs 


tmax 


tmax / 'taver 


10- 4 


3529 sees 


20 


1 


21 


46 


592 sees 


8.4 


10~ 5 


4716 sees 


47 


1 


48 


100 


671 sees 


17.1 


io- y 


7905 sees 


122 


1 


123 


253 


449 sees 


22.0 
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We can see that the problem is unbounded. In Figure [3] the upper image V and the lower 
image T>* for e = 10 -4 are shown. We observe that the vertices of V are almost on a line and 
the lower image T>* is more suitable to illustrate the example. 





Figure 4: The lower image T>* in Example 6.4 for e = 10 , two different view points 



Example 6.4. Now consider d = 11 assets with a given correlation structure and all other 
input parameters as for the first 11 assets in Example 6.3 above. We have N = 2 rf_1 = 
1024 and the number of generating vectors of each cone Kq and i^r(o; n ) for n = 1,...,N 
is d{d — 1) = 110. Consider a basket call option with physical delivery and strike price 
K = ^2i=i(So)i- If at maturity (T = 1 year) the value of the basket of risky assets is greater 
or equal to the strike, i.e., Yll=i(^T)i > K, then one would exercise the option and buy the 



10 risky assets at strike K by delivering K times the bond, i.e. X 



K. 1, 



1) in this 



case. If the value is less, the payoff vector is the zero vector. As the space of eligible assets 
we choose the space spanned by the first three assets, i.e. M = R 3 x {0} 8 . The ordering 
cone is K^ 1 , which is strictly larger than M 3 ^ and generated by 6 vectors. The linear vector 
optimization problem has 3 objectives and a matrix of size 11 272 x 124025, which is sparse 
having 481 288 nonzero entries. For e = 1CP 3 , the computational time of the primal algorithm 
was 1748 seconds. The result is illustrated in Figure [4} As the upper image V is difficult to 
visualize (a polyhedron containing no lines but being "close" to a halfspace) we only provide 
the lower image T>*. 



Example 6.5. Consider cZ = 10 assets with a given correlation structure and all other input 



parameters as in Example 6.4 above. Let X be the payoff of an outperformance option with 
physical delivery as described in Example |6.3| As the space of eligible assets we choose 
the space spanned by the first four assets, i.e. M = M 4 x {0} 6 . The corresponding linear 
vector optimization problem has 4 objectives and a matrix of size 5 126 x 51 300 with 197 638 
nonzero entries. The ordering cone is Kq 1 , which is strictly larger than R 4 and generated 



by 12 vectors. Then, AV@R a (X) with a as in Example 6.4, obtained as the upper image of 
linear vector optimization problem computed with the primal algorithm and e = 10~ 2 , has 
18 vertices and 12 extreme directions. The vertices of V are visualized by a radar chart in 
Figure [5] 



The following numerical example refers to Example 5.2 in the previous section 
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1 




Figure 5: Visualization of the 18 vertices of the upper image V C 
radar chart. 



in Example 6.5 by a 



Example 6.6. Let us consider d = 9 assets with a given correlation structure, all other input 
parameters as in Example 6.4 above (i.e. m = 3), and the same basket option (basket call) 
with payoff X as in Example 6.4 We want to calculate RWC(-X), the relaxed worst case 
risk measure at —X, as described in Example 5.2 with parameter q = 500 for i E {3, ... ,9} 
and zero otherwise. The cone G can be seen as a worst case solvency cone and is chosen 
to be a conservative modification of Kq, where A is replaced by the larger transaction costs 
^wc ~ A -j- 0.2. 





Figure 6: The lower image T>* in Example 6.6 for e 
computed with the dual algorithm. 



10~ 2 (left) and e = 10~ 3 (right) 



RWC(-X) corresponds to an upper good deal bound as it is a relaxed version of the 
set of superhedging portfolios. By considering certain small risks controlled by e and G as 
acceptable, the scalar superhedging price of 34.942464 units of bond is reduced to 34.830995 
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units of bond for RWC{— X). The linear vector optimization problem to calculate RWC{— X) 
has 3 objectives and a matrix B of size 4608 x 36 939, which is sparse with 185 856 nonzero 
entries. The above prices in units of bond were obtained by solving linear vector optimization 
problems with both the primal and dual algorithm for e = 10 -4 . 

The following table shows some computational data and a comparison of the primal and 
dual algorithm. In Figure |6j parts of the results are visualized. 



£ 


variant 


total time 


\S\ 


\T\ 


#LPs 


tmax 


tmax / taver 




primal 


113 sees 


13 


12 


37 


8.8 sees 


3.8 


10" a 


primal 


239 sees 


68 


37 


123 


9.0 sees 


6.8 


10" 4 


primal 


506 sees 


153 


82 


308 


8.8 sees 


8.6 


io- 2 


dual 


86 sees 


7 


20 


37 


5.8 sees 


3.8 


10" 3 


dual 


193 sees 


30 


73 


113 


8.2 sees 


4.8 


10" 4 


dual 


404 sees 


74 


136 


256 


5.8 sees 


5.1 
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